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Summary. The mathematical framework for a-posteriori error estimation of func- 
tionals elucidated by Eriksson et al. [7] and Becker and Rannacher [3] is revisited in a 
space-time context. Using these theories, a hierarchy of exact and approximate error 
representation formulas are presented for use in error estimation and mesh adaptiv- 
ity. Numerical space-time results for simple model problems as well as compressible 
Navier-Stokes flow at Re — 300 over a 2D circular cylinder are then presented to 
demonstrate elements of the error representation theory for time-dependent prob- 
lems. 


1 Introduction 

For better or worse, our physical world is constantly evolving in time. Many impor- 
tant physical phenomena depend fundamentally on time either deterministically or 
through dynamical system behavior. As the complexity of numerical time-dependent 
fluid ffow simulations continues to rapidly increase with advancements in computer 
hardware, the ability to represent, estimate and control numerical errors occurring 
in these time-dependent simulations becomes paramount so that the uncertainty 
and risk associated with using these results in engineering designs becomes accept- 
able. This task is especially difficult for systems such as the Navier-Stokes equations 
where the prospect of controlling pointwise errors is intractable at sufficiently large 
Reynolds numbers. In this latter case, a pointwise error bound exists but the bound 
grows too rapidly with Reynolds number thus rendering the result useless. Even 
so, it still may be feasible to represent, estimate, and control (via adaptivity) the 
error in .certain space and space-time integrated quantities such as body forces, aver- 
ages, and other fluid statistics that can mathematically represented as functionals, 
the mapping of a function space to a single real number. In these special cases, 
the integration process “averages out” the sensitivity of pointwise errors to solution 
perturbations so that useful error bounds are obtained. For a given solution data, 
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u , in m dependent variables, these scalar output functionals will be denoted by 
J(u) : JR m IR. For example, an output functional considered in the next section 
is the space-time average of the i-th solution component in the hypercube centered 
at ( 20 , 2 / 0 , 20 ^ 0 ) 

J[ui) — f f Ui(x — £o, y — 2/o, 2 — 20 , t — to) dx z dt . (1) 

Jto Jdxdxd 

1.1 Computability of Functional Outputs 

To illustrate the benefits of estimating integrated functional outputs, Hoffman and 
Johnson [14] have computed solutions of the backward facing step problem shown in 
Fig. 1 at Re = 2000 using a conforming finite element method with linear elements 
for incompressible flow, see [14] for details. Space-time error estimates were then 
constructed for the functional (1) for various values of the box width d together 
with the time averaging interval to — 9, £i = 10. In velocity and pressure variables 



Fig. 1. Schematic of backward facing step problem showing J(u) averaging box 
position used by Hoffman and Johnson. 


(V,p), the following error estimate for the streamwise velocity component functional 
have been obtained in terms of the backwards in time dual solution (-0, <j>) and the 
element residuals r^,o and rh,i for the (V,p) system 


\J(V,p) - J(V h ,p h )\ < Co\m ||Aor fc ,o(V r kl 'p h )|| +Ci\\D 2 4>\\ \\h 2 r h<0 (V h ,p h )\\ 

+ 0*11011 WhorH^V^W+CsmW 

where ho and h are the temporal and spatial element lengths. Motivated by this 
estimate, stability factors for boxwidths d = 1/8, 1/4, 1/2 have been computed by 
Hoffman and Johnson as summarized in Table 1. 

The growth in stability factor magnitude with decreasing boxwidth size clearly shows 
a deterioration in the error estimate for the functional (1) as a pointwise error 
estimate is approached. In that same work, Hoffman and Johnson also use a drag 
force functional for incompressible flow over a cube geometry attached to a solid wall. 
By using error estimates of the drag functional in mesh adaptation, the mesh is only 
locally refined in a small portion of the domain so that numerical errors in the drag 
force are controlled using a relatively small number of degrees of freedom. For other 
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d 

Ml 

IIWH 

IIWII 

M\ 

1/8 

124.0 

836.0 

138.4 

to 

00 

ih 

1/4 

39.0 

533.4 

48.9 

46.0 

1/2 

10.5 

220.3 

16.1 

25.2 


Table 1. Stability factors for the backward facing step problem for the averaging 
functional (1) for d — 1/8, 1/4, 1/2. 


problems, the overall challenge is not only to find convenient ways to estimate and 
control the numerical error in these functionals, but also to find suitable functionals 
that are both physically meaningful and lead to the resolution of the essential physics 
in the simulation. 

1.2 Estimation of Functional Outputs 

Let u denote the exact solution data for a given problem and Uh a numerical approx- 
imation of this exact solution. The task at hand is the representation, estimation, 
and control of errors in the numerically approximated functional so that given a 
user specified tolerance TOL 


\J{u) - J(u h )\ < TOL 
is achieved using minimal resources. 

At first glance, the inclusion of time derivatives into the formulation appears 
to significantly complicate the construction of error bounds for functionals. Fortu- 
nately, for the class of uniform space-time finite- element methods this is not the case 
since the time dimension is discretized in precisely the same way that any spatial 
dimension is discretized. Specifically, we consider herein a space-time variant of the 
discontinuous Galerkin (DG) method due to Reed and Hill [24] and LeSaint and 
Raviart [19] for hyperbolic problems. In the space-time DG method, the finite el- 
ement discretization in time naturally yields time causality so that discretizations 
can be marched forward in time on a time-slab by time-slab basis for both hyper- 
bolic and parabolic problems. Thus conceptually, the addition of time into the DG 
method does not introduce any new analysis complications. As a practical mat- 
ter, full space-time formulations of DG are quite computationally expensive but 
extremely accurate methods. 


1.3 Preview of Space-Time Error Bounds 

For an introduction to a posteriori error analysis of functionals see the articles by 
Eriksson et al. [7], Becker and Rannacher [3], Giles et al. [9, 10], Johnson et al. [17], 
Prudhomme and Oden [22, 23], Siili [25], the collected NATO lecture notes [1] and 
the multitude of additional references contained therein. In Sect. 2, a brief intro- 
duction to these general theories is given with special attention given to space-time 
generalizations of previous results for the discontinuous Galerkin method previously 
given in Larson and Barth [18] as well as Hartmann and Houston [13]. 
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Our main goal is to present the formulation of error estimates for output func- 
tionals in space-time. To illustrate the general form of these estimates, consider the 
space-time hyperbolic problem with initial data u Q (#) and space- time boundary data 


Cu — s(x,t) m u it + div(/(u)) — s(x,t) — 0, 

in Q x [0, T] 

(2) 

P ~ (v -g(.x,i)) = 0, 

on r x [0, T] 

( 3 ) 

u(x } i ■== 0) uq(x) = 0, 

initial data 

( 4 ) 


where p ± is a characteristic space-time projector such that p~ = 1 at space-time 
inflow and p+ — 0 at space- time outflow. Let JC be a partition of a polygonal do- 
main Q with spatial boundary F into non-overlapping shape regular elements (or 
control volumes) denoted by K , let I n denote an interval in time I n = [t+ , 1 ] , and 

Q n = K x I n a single space-time element with space-time dimension h. Through- 
out, we always consider time integration to a fixed final time T using N discrete 
time slabs such that [0,T] = u£ r “ 0 1 J n . Further, let (a, 6 )q = f Q (a • b) dx dt and 
(a, b)dQ = f dQ {a • b) dxdt. In Sect. 2, the abstract theory is given for space-time 
error representation and estimation of functionals which has historically served as 
the basis for constructing weighted and unweighted estimates of the form (trivially 
extended here to include time): 

• Weighted error estimates for functionals-, Becker and Rannacher [3] 

N 

\J(u)-J(u h )\ < ^^\{r h ,4> ~ ■K h (j>)Qr. + y2(ifc,(/>-7rh<£)eQn| ( 5 ) 

n=0 Q n dQ n 

where tk denotes the residual on element interiors 


r h — Cu h - six, t) 


(6) 


Similarly, jh denotes an interface residual for a space-time element, Q n = Kxl n , 
due to jumps occurring at element interfaces and inflow boundaries, viz. 



6K\r x r 

v - (f(g) ■ n- f(u h ) • n), dKnr x I n 

t n 

, RT, 7i^0 

(uh (t + ) ~ wo (z)), K,n = 0 


( 7 ) 


In the basic error estimate, (j) represents an auxiliary function obtained from 
solving a (backwards in time) dual problem and Hh<t> is any projection of (j) into 
the approximation space containing Uh as discussed in Sect. 2. Suppose J(u) == 
/r?x [o T] 0^ ’ u ) dt ^ or a sufficiently smooth weighting function !^, then the dual 
problem is related to the locally linearized adjoint problem 


together with appropriate boundary conditions. Observe that the functional er- 
ror is a weighted combination of numerical residuals with weighting function 
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• Error estimates for functionals via stability factors, Johnson et al. [17, 14]. The 
computation of the weighted error estimate requires the resolution of a dual 
problem <$> which becomes computationally expensive when performed in space- 
time. In the unweighted estimate of Johnson et al. [17], basic approximation 
theory results are applied to the dual problem appearing in (5) thus yielding 

( N~i 

y. + hS l/2 ib^iidQ n ) 2 

n— 0 Q n 

for 1 < s < p -f 1. The dual problem is then replaced with a stability constant 
satisfying 

||<£||h s (x?x[o,t]) < Cstab 

thereby producing the computable estimation formula 

( N - 1 

E + hS ~ 1/2 \\MM 2 

n—0 Q n 

The basis idea is that it may be feasible to catalog stability constants C sta b for 
various flow problems so that computation of the dual problem can be avoided. 
The simplifying assumptions utilized here lead to an unavoidable loss in sharp- 
ness when compared to the weighted error estimate. 

Next, we briefly review the relevant theory for representing the space-time error 
in functionals which serves as the primogenitor for both of the previously given 
a-posteriori estimates. 



2 Error Estimation of Functionals 

Consider once again the model space-time problem (4) in a domain Q C IR d with 
boundary T. In this domain, let V® )33 denote a mesh dependent broken space of 
piecewise polynomials of complete degree p in each Q n with no continuity between 
space-time elements 

Vl P = {w : w\ Q € Vp(Q) , VQ e K x U^ 1 J n } . (8) 

To illustrate the abstract framework for error estimation of functionals, we consider 
the discontinuous Galerkin (DG) finite element method introduced by Reed and 
Hill [24] and LeSaint and Raviart [19] as analyzed by Johnson and Pitkaranta [16] 
and further refined for nonlinear conservation laws by Cockburn et al. [4, 5]. In 
presenting the abstract theory, one quickly sees that the techniques used here can 
be easily applied to any method that can be written in weighted residual form. 

Space-Time Discontinuous Galerkin FEM. Find Uh £ Vh lP such that 


where 


&DG(u h , w) — F(w ), Vw e Vh,p 


( 9 ) 
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Bd G( v>h,w) - F(w) = y; j y A I (-w.t-Ufc . fi(u h ))dxdt 

n=0 \q« VQn 


+ 

+ 

+ 


/ / u>- • h(n\ (uh)-, (uh)+) dsdt 

Ji n J (dKndQ n )\r 

/ / u;_ ♦ h(n; (u&) _ , g) ds dt ] 

i(OKnaQ^)nr ' 

[ (w(t" +1 ) Wh (C +1 ) - t»(4) u„(t*)) dx 

f (u;(£l) wa(^) - w(i+) uo)) ch 

J f2,n^0 


where /i(n; < u_, / iq-) is a numerical flux function such that h{n\ u, u) — ni f%(u) and 
h{n\ u- t u+) == — fe(— -n;w+,u^). 


2.1 DG FEM Error Representation: The Linear Case 

The objective is to estimate the error in a user specified functional J(u), the mapping 
of a function space to a single real number. In this work, we consider functionals 
that can be expressed as a weighted integration over the space-time domain 

J^(u) = f [ ip - N(u) dxdt 

Jo Ja 

or a weighted integration on the boundary of the space-time domain 

J^(u)~ ip • N(u) ds 

Jd(S 2 x[ 0 t T\) 

for some user specified weighting function : IR^ 3R m and linear /nonlinear 
function N(u) : IR m i — > lR m . By an appropriate choice of ip(x) and N(u), it is 
possible to devise functionals of practical engineering use, e.g. lift and drag forces 
on a body, stress intensity factors, average quantities, etc. 

Let J5 dg(*, 0 denote a bilinear form for the discontinuous Galerkin method and 
J(-) a linear functional. In the following derivations, Wh denotes any suitable projec- 
tion operator (e.g. interpolation, L 2 projection) into yt . Begin by introducing the 
primal numerical method assuming all boundary conditions are weakly enforced. 

Primal numerical problem: Find Uh € V^ )P such that 

BT>G(uh, w) = F(w) V w € Vl P 

with the Galerkin orthogonality condition 

B&g{u - u h , w) = 0 V w € Vh tP . 

Next, we introduce the auxiliary dual problem utilizing infinite-dimensional trial 
and test spaces. 

Dual problem: Find $ € V B such that 

£dg(w,£) == J(w) \f w £V B . 
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An exact error representation formula for a given functional J(-) results from the 
following steps 

J(u) - J{u h ) — J (u - u h ) 

— Br>o(u — Uh , $) 

== Bv>g(u ^ Uh , # — Kh&) 

= I?Dg(u, # — TT/i^) — #DG (Uh, ^ — 7Th^) 

— F($ ■— 7T^#) — J5dg(^/i> — 7Th^) 

so in summary 

J(u) — J(uh) = F($ — 7Th$) — 5 dg(^, $ — 

Notably absent from the right-hand side of this equation is 
exact solution u. 

2.2 DG FEM Error Representation: The Nonlinear Case 

Let Bdg(-> ') denote a semilinear form for the discontinuous Galerkin method and 
J(-) now a nonlinear functional. To cope with nonlinearity, mean- value linearization 
is employed 

£>dg(u, w) = Bdg (Uhi uj) -f Bdg (uh, u; u — Uh, u>) V ty G V B (11) 

J(u) = J(uh) + T(uh, u; u — Uh) . (12) 

For example, if B(u, w) == (£u,u;) for some nonlinear differential operator £ then 
for u> G V 

S(u,u;) = £(uh,u>) + ^ J C, u {u(9))d0 (u — Uh),u^ 

— #(uh, u;) + (£ >u (u - Uh) , yj) 

— B(uh , u;) -f B(uh, u; u - Uh, vi). 

with u(0) = Uh + (u — Uh) For brevity, the dependence of £> on the path integration 
involving the exact solution u will be notationally suppressed. We then proceed in 
the same fashion as in the previous example. 

Primal numerical problem: Find Uh € V B P such that 

Bdg{uh,w) = F(w) V w G Vh lP (13) 

with orthogonality condition for the linearized form 

Bvg(u — Uh, u;) 0 V w € V B P . 

A mean-value linearized dual problem is then introduced which utilizes infinite- 
dimensional trial and test spaces. 

Linearized dual problem: Find $ G V B such that 

Bvq{W)<F) — J(w) V w G V B . 


(linearity of M) 

(dual problem) 
(orthogonality) 
(linearity of B) 
(primal problem) 

7Th$) . (10) 

any dependence on the 


(14) 
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An exact error representation formula for a given nonlinear functional J(-) then 
results from the following steps 


J (u) - J (u h ) = J{u - u h ) 

— Bdg(u — Uhy$) 

— Bdg(u - Uh,$ — 7V h 0) 

— 0 — lTh&) “ BdG (*%, $ ~ TTh,#) 

= F{0 — Kh&) — BuG{uh) 0 -~ 


(mean- value J) 
(dual problem) 
(orthogonality) 
(mean- value B) 
(primal problem) 


so in summary 

J(u) - J(u h ) = P(# - **#) - &a(w h , # - 7T^) . (15) 

Note that although Eqns. (10) and (15) appear identical, mean-value linearization 
introduces a subtle right-hand side dependency on the exact solution in Eqn. (15). 


3 Computable Error Estimates and Adaptivity 

Computationally, the error representation formulas (10) and (15) are not suitable 
for obtaining computable a posteriori error estimates and use in mesh adaptation. 

• The function 0 — i Th0 is unknown where 0 € V B is a solution of the infinite- 
dimensional dual problem. 

• The mean- value linearization used in the linearized dual problems (14) requires 
knowledge of the exact solution u. 

Various strategies which address the numerical approximation of 0 are discussed 
in Barth and Larson [2], e.g. postprocessing, higher order solves, etc. Due to Galerkin 
orthogonality, the dual problem in the discontinuous Galerkin finite element method 
must- be approximated in a larger space of functions than that- utilized in the primal 
numerical problem. For purposes of the present study, this is achieved in the discon- 
tinuous Galerkin method by solving the dual problem using a polynomial space that 
is one polynomial degree higher than the primal numerical problem, viz. if Uh € 
then g V® p+1 . 

In the present study, the mean-value linearization depending on the states u and 
Uh is replaced by the simpler Jacobian linearization evaluated at the numerical state 
Uh . This is not the only practical choice. In Barth and Larson [2], a more sophisti- 
cated technique involving the postprocessing of primal data and the approximation 
of the mean-value linearization by numerical quadrature is employed in computa- 
tions. 

3.1 Direct Estimates 

Whin written in global abstract form, the error representation formula does not 
indicate which elements in the mesh should be refined to reduce the measured error 
in a functional. By applying a sequence of direct estimates, error bounds suitable 
for adaptive meshing are easily obtained. The goal in constructing these estimates 
is to estimate the local contribution of each element in the mesh to the functional 
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error. This local cell contribution will then be used as an error indicator for choosing 
which elements to refine or coarsen in the adaptive mesh procedure. 

| J(u)— J(uh)\ = \Bi)G(uh , £ ~ Kh$) — F(# - | (error representation) 

N - 1 

= | i^ DG >Q n ~ Kh$) — Fq n (# — Kh&))\ (element assembly) 

n - 0 Q n 
N-X 

< r T: 1 (u h ,$- 7T h $) - Fq » (# - 7^#)) | (triangle inequality) 

77.— 0 Q™ 

(16) 


where •) and Fq« (•) are restrictions of Bdg(-, *) and FQ to the partition 

element Q n . The basic definition of the discontinuous Galerkin semilinear form given 
in (10) shows one possible element assembly form but this is not a unique repre- 
sentation. For example strong and weak forms of the semilinear operator j^dg^*,-) 
yield differing assembly representations. For the discontinuous Galerkin method, the 
error representation formula together with (10) for a single element Q n yields 


£dg,Q" Wh#) - F Q -n (# - 7 Th$) = - / 7 u h • ($ - 7 Th$),t dx dt 

Jin Jk 

- [ I fi{uh) ■ (# - dx dt 

dj* d/c 


+ 

T 

+ 


/ / (<£-7r^)_ * h(n;(n/>)- 5 CuQ-f )dsdt 

d/^ JdK\r 

/ / (#-717^)- -/i(n;(n4- } p)dsdt 

d/n daicnr 

L 
L 
L 


($~7T h $)(tt +1 )‘U{tV- l )d X 


K 


i<,n#0 


K,n-0 


(#*-7 Thtf)(4)-«(C)dp 

3 

(# - 7T^)(t+) * n 0 ) dx 


or a weighted residual (strong) form can be obtained upon integration by parts 
BDG t Q n ('Uh)& — TTh,#) — Fq^ — iijiF) == / {<£ — 7rad^) * dx dt 

Jq n 

~f / / (# - 7r h ^)~ • (M n ; (^)~ } (^)+) - /( n ; (^)-)dsdt 

dp* JdK\r 

+ / (^-7r/ l #)_-(/i(n;(nh)- 5 p)-/(n;(n h )„)dsdt) 

dpv JdKnr J 

+ [ ($-1 r^)(4) • [u]*+ dx 

JK t n^0 

+ [ (#- 7r h #)(t° )-(uh(t+)~ u 0 )dx (18) 

JK,n=0 


This latter weighted residualform #Q n (*> *) ~ Fq^Q is preferred in the error 
estimates (16) since the individual terms represent residual components that vanish 
individually when the exact solution is inserted into the variational form and a 
slightly sharper approximation is obtained after application of the triangle inequality 
in (16). 
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3.2 Adaptive Meshing 

Motivated by the direct estimates (16), we define for each partition element Q n the 
space-time element quantity 

rjQn = - TTfe#) “ Fq™ (# — 7T/,.#) , (19) 

so that \r)Qn | serves as an adaptation element indicator such that 

|J(«)-JK)|< (20) 

7i=0 Q™ 

and the sum of rjQn over space-time provides an accurate adaptation stopping criteria 

ij(«)-j(« h )i = iy;X>"i - (2i) 

n=0 Q n 

These quantities suggest a simple mesh adaptation strategy in common use with 
other indicator functions: 

Mesh Adaptation Algorithm 

(1) Construct an initial space-time mesh. 

(2) Compute a numerical approximation of the primal problem on the current space- 
time mesh. 

(3) Compute a suitable numerical approximation of the infinite-dimensional dual 
problem on the current mesh. 

(4) Compute error indicators, rjqn . 

(5) if( I £n=o I < tol ) stop 

(6) Otherwise, refine and coarsen a specified fraction of the total number of elements 
according to the size of |??|q, generate a new mesh and GOTO 2 

3.3 Numerical Example: Error Estimates for Time-Dependent 
Hyperbolic Problems 

Numerical error estimation results for stationary hyperbolic problems using the dis- 
continuous Galerkin method and the error estimates described above are given by 
the present author in [18, 2] and later revisited by Hartman and Houston in [13] . 
To demonstrate the application of the error representation formula (10) and the 
error bound (16) to time-dependent hyperbolic problems, discontinuous Galerkin 
solutions have been obtained for the following 2-D advection problem in the square 
domain Q G [—1, l] 2 


u )t + A • Vu = 0 
u(x } t) — 0 

u(x, 0) — uq(x) 


in Q x [0, T] , 
on dQ x [0, T] 
initial data 


( 22 ) 


with circular advection field A = (— y 1 x) T i final time T — 1.15, and C°° initial data 
uo(x) = cr(2/10; y/Jx -Xo) 2 + (y- yo) 2 ) , (xo,yo) = (7/10, 0) 
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with 

v{ro;r) = j el+r g ^ ^ ■ ( 23 ) 

An output functional was also devised consisting of a C°° space-time averaging ball 
of radius 1/4 located at (x — (.30, .26, .64), i.e. 

J{u) — f f (j(l/3; y/(t — tf) 2 + (x — Xf) 2 + (y - y f ) 2 ) dxdt . 

Jo Jn 



Fig. 2. Primal solution of the circular advection problem at t — .45. Shown are 
a carpet plot of the primal solution using linear space-time elements (left) and a 
graph of the accumulated functional error, J(u) — J(uh), during the backwards in 
time dual solution integration (right) . 



Fig. 3. Dual solution of the circular advection problem at i — .45. Shown axe a 
carpet plot of the dual solution 4>h using quadratic space-time elements (left) and 
a carpet plot of the dual solution defect <j)h — ?r h<j>h (right) with tth a projection to 
space-time linear functions. 


The numerical experiment then consists for the following steps 
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1. Numerical solution, Uh, of the primal (forwards in time) problem using the 
discontinuous Galerkin method (9) with linear space-time elements, see Fig. 
2(left). 

2. Numerical solution, fa, of the dual (backwards in time) problem using the 
discontinuous Galerkin method (9) with quadratic space-time elements, see Fig. 

3. 

3. Approximation of the infinite-dimensional dual problem by the quadratic space- 
time element approximation, <j> « fa. 

4. Accumulation of the functional error, | J(u) — J{uj) |, using the error represen- 
tation formula (10) and the error bound formula (16), see Fig. 2 (right). 

Note that although we plot in Fig. 2 the accumulation of terms appearing in the 
error estimation formulas (10) and (16) during the backwards in time integration, it 
is only the final termination value at time zero that provides a bound for the desired 
functional error. As expected, the error representation formula very accurately re- 
produces the exact functional error for this linear differential equation since the only 
source of approximation in this formula is the replacement of <f> by a higher order 
quadratic space- time element approximation fa. The graph in Fig. 2 corresponding 
to Eqn. (16) lies strictly above the error representation formula since the absolute 
value occurring in the Eqn. (16) precludes any interelement cancellation of error 
terms. In this particular example, this lack of interelement cancellation in Eqn. (16) 
also gives rise to a different slope and a continuous deterioration in Eqn. 16 with 
increasing functional measurement time. 


4 Extension to the Time-Dependent Navier-Stokes 
Equations 

Next, we consider the error estimation of functionals applied to the time-dependent 
compressible Navier-Stokes equations for a perfect gas with fluid density, velocity, 
pressure, and temperature denoted by p, V y p and T 

U,t + - F™ in n x [0, T\ (24) 

with implied summation on repeated indices, 1 , . . . , d, and 



( p \ 


pVi \ 


0 \ 



pV i 


pViV\ + Sup 

t-ivis 1 

’ Fi ~ We 

n< 


u — 

pV d 

, 4 nv = 

• 

pViVd + SidP 

TZi 

(25) 



\pE ) 


V ViiE + p) J 


\TijVj + 



T = £ (W + ( VK ) T - f( V • V W ( 26 ) 

where 7, Re, and Pr denote the non-dimensional ratio of specific heats, Reynolds 
number and Prandtl numbers, respectively. 

The present implementation of the discontinuous Galerkin method utilizes a 
change of dependent variables, u u, where v are the so-called symmetrization 
variables for the Navier-Stokes equations as described in [15]. Using this change of 
variables, viscous fluxes are written in chainrule form 



Space-Time Error Representation and Estimation in CFD 


13 


F?' ls - Mij v tXj , Mij € IR mxm . 

The extension of the discontinuous Galerkin method to the time-dependent Navier- 
Stokes equations follows the original symmetric interior penalty (SIP) work of Dou- 
glas and Dupont [6] for elliptic and parabolic problems. The SIP formulation has 
been recently applied to the steady Navier-Stokes equations in [11, 12]. The present 
SIP formulation and subsequent error estimation calculations accommodates the 
full time-dependent compressible Navier-Stokes equations. For brevity, the method 
is stated here in strong form. In practice, the method should always be implemented 
in weak integrated-by-parts form so that exact discrete conservation is assured even 
in the presence of inexact numerical quadrature. 

Space-Time Discontinuous Galerkin SIP Navier-Stokes FEM. Find vn € such 
that 

N^l 

Y Y ^DG,g^ ( vh , w) = 0 , Vu> € Vh iP (27) 


Bvo,q~(v,w) = J J w ■ (u,. + Fi™ - F & ) 


dtdx 


+ 

+ 

+ 


+ 


+ 

+ 


f f ui(x~) • (h(n; u_) — rii F* nv (u_)) dt dx 

JdK\r Ji n 

f [ vj(x-)- m ( F ? v wal1 - Ft")(v - )) dt dx 

JdKnr w&n J i* 

f f w(x~) • (/i(n; goo^vJ) — Ui i^ nv (u_)) dtdx 

JdKnr {&r{ield Jin 

[ f 7; w ( x -) * [ n i J p i ls ]xt dtdx 
JdK\rJin 2 

[ [ w(x-) • rii (Fi is (g N ) - FP s (v _)) dtdS 

JdKnr N Jin 

/ / o * n i Mij ( x ~ ) W ,Xj ( x - ) dt dx 

JdK\r Jin 2 

/ / (9d — v (x~)) • niMij(x-)w iXj (x~) dtdx 

JdKnr D Ji n 

f f (Kp 2 /h)xt w(x-)) • rii rij Mij [v]®* dt dx 
JdK\r Jin 

/ / (Kp 2 /h) x _ w(x _)) • rii rij Mij (gn — v(x~))) dt dx 

JdKnrv Jin 

L 


K,n^ 0 


u>(t+) * [u(u)] £ i dx 


f W (J+) ' (u(u(£° J) — uo) dx 
JK,n = 0 


where k is an (9(1) parameter, {•)£* = ((•)»_ + (')x + )/2 and g oo, 9n and go denote 
far field, Neumann, and Dirichlet data, respectively. This formulation immediately 
gives a convenient error representation formula of the same form as (15) 


N-l 

J(u(vh)) “ J(u(v)) = Y, 

n - 0 K 


( 28 ) 
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and an error estimate suitable for mesh adaptivity of the same form as (16) 

N - 1 

\J(u(v h )) “ J(u{v)) j < ^ ^ |BdG,Q^ (Vh> <t> - 7Ch<f>) :| • (29) 

n=0 K 

4.1 Numerical Example: Time-Dependent Compressible 
Navier-Stokes 2D Cylinder Flow at Re — 300 

The accurate prediction of lift and drag forces on a cylinder undergoing vor- 
tex shedding remains an important and difficult test problem in numerical time- 
dependent Navier-Stokes simulations [21, 8, 20]. In the present study, subsonic flow 
at Mach=l/10 and Re=300 has been computed over a 2D cylinder for a length of 
time sufficient to establish periodic vortex shedding. A numerical solution suitable 
for assessing error estimates and adaptation was computed using 12 K linear space- 
time elements and a time step size resulting in approximately 48 time steps per 
oscillation period. Another extremely accurate “reference” solution was also com- 
puted using 40 K quadratic space-time elements and 50% reduced timestep, see Fig. 
4. 



Fig. 4. Navier-Stokes solution at the non-dimensional time t — 635 computed on the 
reference 40K element mesh using Vi space-time elements. Presented here are veloc- 
ity contours (left) and logarithmically scaled vorticity magnitude contours (right). 


A measurement functional was then devised consisting of the the non-dimensional 
drag force integrated over the cylinder surface and averaged over four drag oscil- 
lation periods, t € [605,700], determined from the reference solution with time 
non-dimensionalized here using freestream sound speed and cylinder diameter. 

J(u) — y- 2 7 ~T f ( Fdrag(u) dxdt . (3<5) 

l/^poo Voo tb — to Jt a J cyl wall 

Following the same procedure as outline earlier in Sect. 3.3, a dual solution corre- 
sponding to the time-averaged drag functional was obtained on the 12 K mesh us- 
ing quadratic space-time elements. This solution serves as an approximation of the 
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infinite-dimension dual solution <f>. Figure 5 provides a time snapshot of the quadratic 
space-time element dual solution, <f>h , and the defect dual solution, <j>h — ^k<j>h , where 
7T/t is the 1/2 space-time projection from quadratic elements to linear elements. A 



Fig. 5. Locally linearized dual Navier-Stokes solution corresponding to the drag 
coefficient functional with local linearization about Vi primal space-time data on a 
12K mesh. Presented here are dual solution contours at the same approximate non- 
dimensional time as Fig. 4 showing the ^-momentum component of the dual solution 
4> (left) and the x-momentum component of the dual solution defect <j> — 7 Th<t> (right). 


rather striking feature of the dual solution as is illustrated in Fig. 5 is the local- 
ized support of the defect <f> — 717 ^ when compared to <j> itself. Since only the defect 
<l> — 7Th(j> appears in the error representation formula (28), the resulting mesh adaptiv- 
ity becomes very localized as well. Although it is possible to adapt the mesh during 
every time step of the calculation based on error indicators derived from (29), in 
the present calculations the mesh indicators have been averaged over four oscillation 
periods so that a single (static) mesh adaptation is performed, see Fig. 6. Figure 7 




Fig. 6. Closeup of the original 12 K element mesh (left) and refined 15 K element 
mesh (right) obtained from error indicators averaged over four oscillation periods. 
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shows the lift and drag coefficient histories on all meshes. A comparison of drag co- 
efficient histories shown in this figure reveals a dramatic improvement in oscillation 
amplitude and mean value with only modest mesh adaptivity. Table 2 quantifies the 




Fig, 7. Lift and drag coefficient history for 2D Navier-Stokes cylinder flow at 
Re— 300 and Mach=0.1 using 12 K V\ space-time elements, 15 K Vi adapted space- 
time elements and 40 AT V 2 space-time elements. Shown are the lift and drag coeffi- 
cient histories (left) and an extreme closeup of the drag coefficient histories in the 
functional interval (right). 


accuracy of the computed functional error using the 40 K element reference solution 
as the “exact” solution together with the estimated values from (28) and (29). The 


^elements 

1 J(u) - J(u h ) I 

Eqn. (28) 

Eqn. (29) 

12K 

15K (adapted) 

.0107 

.0020 

.0183 

.0024 

.1098 

.0139 


Table 2. Error estimates for the 2D cylinder on meshes with 12 K and 15K - adapted 
elements. 


tabulated results show some lack of sharpness in the computed error representation 
formula (28) that would be further amplified in the estimate (29). The improvement 
with adaptive mesh refinement may indicate that this source of error originates from 
linearization error resulting from the use of the Jacobian linearization as a replace- 
ment for the mean- value linearization in (11) and (12). Another source of error may 
originate from the specification of the four period averaging interval in the functional 
definition. These issues will be addressed in future expanded work for this problem. 


5 Concluding Remarks 

As one might expect, the extension of a-posteriori error estimation to space-time 
is almost transparent for finite element methods that permit discretization in time 






Space-Time Error Representation and Estimation in CFD 


17 


using the same finite element method. Even so, the implied added computation is 
significant since for time-dependent problems a backwards in time dual problem 
must be solved or approximated. Many open problems exist in this problem area: 

• Efficient techniques for approximating the dual in space-time problem. 

• Computable functionals for time-dependent problems that do not lead to dual 
problems that become large or unbounded in reverse time. 

• Techniques for estimating the error in general L p norms. 
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